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Abstract 

Self-sustained musical instruments (bowed string, woodwind and brass instruments) can be mod- 
eled by nonlinear dynamical systems. Among these instruments, flutes and flue organ pipes present 
the particularity to be modeled as a delay dynamical system. In this paper, such a system, a 
toy model of flute-like instruments, is studied using numerical continuation. Equilibrium and pe- 
riodic solutions are explored with respect to the blowing pressure, with focus on amplitude and 
frequency evolutions along the different solution branches, as well as "jumps" between periodic so- 
lution branches. The influence of a second model parameter (namely the inharmonicity) on the 
behaviour of the system is addressed. It is shown that harmonicity plays a key role in the presence 
of hysteresis or quasi-periodic regime. Throughout the paper, experimental results on a real instru- 
ment are presented to illustrate various phenomena, and allow some qualitative comparisons with 
numerical results. 



1 Introduction 

Sound production in self-sustained musical instru- 
ments, like bowed string instruments, reed instru- 
ments or flutes, involves the conversion of a quasi- 
static energy source (provided by the musician) 
into acoustic energy. This generation of auto- 
oscillations necessarily implies nonlinear mecha- 
nisms. Thus, even the simplest models of self- 
sustained musical instruments should include non- 
linear terms |27] , 

We focus in this paper on flute-like instruments. 
Many works with both experimental and model- 
ing approaches have highlighted the wide variety 
of their oscillation regimes, and aim to explore the 
complexity of their dynamics |201 [2T] 123] , How- 
ever, to the authors knowledge, none study has 
ever investigated the determination of solution 
branches using numerical continuation, and the 
analysis of jumps between branches. This is done 
in this paper. 

A numerical continuation approach has recently 
proved to be useful in the context of musical in- 



struments, especially for reed instruments (like 
the clarinet) [3]. However, unlike reed instru- 
ments, flute-like instruments are delay dynam- 
ical systems (see section [2]), which complicates 
the model analysis, and prevents the use of the 
same numerical tools (as, for example, Auto |18| 
or Manlab \19\). 



1 Corresponding author: terrien@lma.cnrs-mrs.fr 



Using a numerical continuation software ded- 
icated to delay dynamical systems, we study in 
this paper a dynamical system inspired by the 
physics of flute- like instruments. We call it a 
toy model, since it is known in musical acoustics 
that more accurate (yet more complex) models 
should be considered. We investigate throughout 
the paper the diversity of oscillation regimes of 
the toy model, as well as bifurcations and jumps 
between branches. We particularly focus on the 
nature of the solutions (static or oscillating), and 
in the later case, on the evolution of frequency 
and amplitude along periodic solution branches. 
Indeed, in a musical context, periodic solutions 
correspond to notes, and oscillating frequency and 
amplitude are respectively related to the pitch 
and the intensity of the note played. Compar- 
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isons between the obtained bifurcation diagrams, 
experimental data and time-domain simulations 
highlight the valuable contribution of numerical 
continuation. Moreover, surprisingly enough, the 
qualitative behaviour of the toy model displays 
close similarities to that of a real instrument. 

The paper is structured as follows. The dynam- 
ical system under study is presented in section [2j 
In section O we present the results of numerical 
continuation of the branches of static and peri- 
odic solutions and their stability analysis, with 
focus on amplitude and frequency evolutions. The 
study of "jumps" between different periodic solu- 
tions branches is presented in section HI Finally, 
we discuss, in section the influence of a parame- 
ter related to the instrument makers' work on the 
characteristics of the different oscillating regimes. 
Qualitative comparisons with the sound produced 
by real instruments and with the experience of an 
instrument maker are presented. 

2 System studied (toy model) 

In this paper, we study the following delay dy- 
namical system: 

z(t) = v(t - t) 

p(t) = a -tanh[z(t)} (1) 
V(oj) = Y(u) ■ P(u) 

where lowercase variables are written in the time 
domain, and uppercase variables in the frequency 
domain. The unknowns are z, v, and p, while a 
and t are scalar parameters. Y is a known func- 
tion of the frequency, and will be detailed later. 
We first briefly describe the mechanism of sound 
production in flute-like instruments, and then we 
precise to which extent this system can be inter- 
preted as an extremely simplified model of this 
kind of musical instruments. 

2.1 Sound production in flute-like in- 
struments 

Since the work of Helmholtz [T], a classical ap- 
proach consists in modeling a musical instrument 
by a nonlinear coupling between an exciter and a 
resonator. In other wind instruments, the exciter 
involves the vibration of a solid element: a cane 
reed in the clarinet or the saxophone, the musi- 
cian's lips in the trumpet or the trombone... In 
flute-like instruments (which includes recorders, 
transverse flutes, flue organ pipes...), the exciter 
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Figure 1: Recorder section, and simplified repre- 
sentation of the jet oscillation on both sides of the 
labium. 

consists in the oscillation of the blown air jet 
around an edge called labium (see figure [T]) [TT] . 

The jet-labium interaction constitutes the 
pressure source which excites the resonator 
formed by the pipe, and thus creates an acoustic 
field in the instrument. In turn, the acoustic field 
disturbs the jet at the channel exit. As the jet is 
naturally unstable, this perturbation is convected 
and amplified along the jet, from the channel exit 
to the labium, which sustains the oscillations of 
the jet around the labium, and thus the sound 
production [2|. 

The convection time of perturbations along the 
jet introduces a delay in the system, whose value 
is related to the convection velocity c v of these 
hydrodynamic perturbations. Theoretical [5l [12] 
and experimental [3] studies have shown that c v 
is related to the jet velocity Uj (itself related to 
the fact that the musician blows hard or not in 
the instrument) through: 

c v ^0.5-Uj. (2) 

Noting W the length of the jet (i.e. the distance 
between the channel exit and the labium, high- 
lighted in figure [TJ , an approximation of the delay 
value is given by: 

= W„ W 
T c v ~ 0.5 -Uj' U 

2.2 System studied: a toy model of 
flute-like instruments 

This mechanism of sound production can be mod- 
eled by a nonlinear oscillator, such as the one rep- 
resented in figure [2] [TT| . This very basic modeling 
includes the three following elements: 
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Figure 2: Basic modeling of sound production 
mechanism in flute- like instruments, as a system 
with a feedback loop [2J. 

• a delay, related to the convection time of the 
hydrodynamic perturbations along the jet. 



a nonlinear function, 
labium interaction. 



representing the jet- 



• a linear transfer function representing the 
acoustic behaviour of the resonator. 

Therefore, system ([I]) can be seen as a very basic 
model of flute-like instruments, and thus can be 
related to the nonlinear oscillator represented 
in figure [2j Variables z, p and v are then 
respectively related to the transversal deflection 
of the jet at the labium, to the pressure source 
created by the jet-labium interaction, and to the 
acoustic wave velocity at the input of the pipe 
(see figure [T]) . The scalar parameters a and r are 
respectively associated to the amplification of the 
perturbations along the jet, and to the convection 
delay defined by equation ([3]). 

Using a modal decomposition of the resonator 
acoustical response, the transfer function is writ- 
ten in the frequency domain as a sum of m modes 
(as it is done, for example, in the work of Silva 

my. 
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(4) 



where a^, and are respectively the modal 
amplitude, the resonance pulsation and the 
quality factor of the k th mode. 

As it does not take into account the following el- 
ements, this representation is very simplified com- 
pared to the most recent models of flute-like in- 
struments |2|: 



a precise modeling of the aeroacoustic source, 
which includes, in the latest models, a time 
derivative of the delayed term 

nonlinear losses, due to vortex shedding at 
the labium 1221. 



• an accurate description of the resonator: in 
theory, the modal decomposition in equation 
(|4|) includes an infinite number of modes. For 
sake of simplicity, we only retain in this paper 
at most the first two modes of the pipe. 

2.3 Rewritting as a first-order system 

To be analysed with classical numerical continua- 
tion methods for delay differential equations |14j . 
system ([T]) should be rewritten as a classical first- 
order delay system x = f(x(t),x(t — r),7). In 
such a formulation, x is the vector of state vari- 
ables and 7 the set of parameters. To improve 
numerical conditioning of the problem, a dimen- 
sionless time variable and a dimensionless convec- 
tion delay are defined : 



t -- 

T 



OJ\t 
- UJ\T 



(5) 



with oj\ the first resonance pulsation (see equation 
Rewritting system ([T]) finally leads to the 
following system of 2m equations, where m is the 
number of modes in the transfer function Y(u) 
(see appendix lAl for more details). 
V i = [1, m] (where i is an integer): 



v(t) 






k=l 


yit) 


m 

= X>(*) 
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Vi{t) 


= y l {t) 


m{t) 





(6) 



{l-tanh 2 [v(t-f)]} 



Vi{t) 



yi(t). 



System ([6]) is studied throughout this paper. Nu- 
merical results are computed using the software 
DDE Biftool [6], which performs numerical con- 
tinuation of delay differential equations using a 
prediction/correction approach |14] , Periodic so- 
lutions are computed using orthogonal collocation 
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|24| . Typically, we use about 100 points per pe- 
riod, and the degree of the Lagrange interpolation 
polynomial is 3 or 4. 

3 Periodic regimes emerging 
from the equilibrium solution 

3.1 Branch of equilibrium solutions 

In flute-like instruments, the delay r is related to 
the jet velocity (see equation ©), and therefore 
to the static pressure into the musician's mouth 
through the stationary Bernoulli equation. To an- 
alyze the system ([6]), it is thus particularly inter- 
esting to choose this variable as the continuation 
parameter. 

Regardless of the parameter values, it is obvious 
that system ([6]) has a unique equilibrium solution 
defined by: 

C Vi= [l,2,..,m] : 

{ v l (t) = (7) 

{ yS) = o 

A stability analysis of this static solution is per- 
formed both through computation of the eigenval- 
ues of the linearized problem and the analysis of 
the open- loop gain (see appendix IB1 for more de- 
tails). Moreover, this linear analysis allows to dis- 
tinguish two kinds of emerging periodic solutions: 

• those resulting from the coupling between an 
acoustic mode of the resonator (a particular 
term in the sum defining the transfer func- 
tion Y(lo) in equation (j3J)) and the first hy- 
drodynamic mode of the jet (corresponding 
to n = in equation (124p of appendix IB"]) . 

• those resulting from the coupling between an 
acoustic mode of the resonator and an higher- 
ranked hydrodynamic mode of the jet (corre- 
sponding to n > in equation (|24p ). 

Indeed, for flute-like instruments, an auto- 
oscillation results from the coupling between 
an acoustic mode of the resonator and an hy- 
drodynamic mode of the jet [2]. The first case 
n = corresponds to the "standard" regimes 
in flute- like instruments. On the contrary, the 
second case n > corresponds to the so-called 
"aeolian" regimes. In the context of musical 
instruments, it corresponds to particular sounds, 
like for example those produced when the wind 
machine of an organ is turned off leaving a key 



pressed. 

We focus here on these two kinds of periodic 
solutions, and we study particularly the amplitude 
and frequency evolutions along the branches. 

3.2 Branches of periodic solutions: ae- 
olian and standard regimes 

For sake of clarity, we consider in this section 
a transfer function containing a single acoustic 
mode (i.e. m = 1 in the expression of Y(u), de- 
fined by equation @). The addition of a second 
acoustic mode will be discussed in the next sec- 
tion. 

Numerical continuation of the different periodic 
solution branches allows to display the bifurcation 
diagram showed in figure [3j representing the am- 
plitude of the oscillating variable v(t) defined in 
system ©, with respect to the dimensionless de- 
lay f. It is useful to keep in mind that blowing 
harder makes the jet velocity Uj increase and thus 
t decrease (according to equation ©). 

The stability of each branch is adressed by com- 
puting the Floquet multipliers: since all of them 
lie within the unit circle, it can be concluded that 
all the branches are stable |10j . 

As shown in figure [3j the use of equation (|24p 
(appendix [B]) allows to distinguish the different 
kind of periodic regimes, and highlights that the 
only standard regime (related to the first hydro- 
dynamic mode of the jet) is located in the part of 
the graph where oj\T is lower than r* = 2. The 
other branches correspond to aeolian regimes re- 
lated to the second and the third hydrodynamic 
mode of the jet (respectively n = 1 and n = 2). 

A larger range in f in figure [3] would reveal 
other branches of aeolian sounds. They corre- 
spond to the infinite series of aeolian instabilities 
highlighted for each acoustic mode of the transfer 
function Y(uj) in appendix [Bj 

3.2.1 Amplitude evolutions 

To study flute-like instruments, it is useful to de- 
fine the dimensionless jet velocity 9 [7]: 



where W is the jet length (see figure [T]), and / the 
oscillation frequency. 
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Figure 3: Bifurcation diagram of system (J6j. Ab- 
scissa: dimensionless delay oj\t. Ordinate: am- 
plitude of the oscillating variable v(t). Symbols 
x represent the stable parts of the branches. Pa- 
rameters value: m = 1, a = 10 ; a\ = 70 ; 
uj\ = 2260 ; Qi = 50. Values of n are computed 
using equation (j24"l) in appendix [Bj 

Indeed, the representation, in ngure[U of the os- 
cillation amplitude of the variable v(t) as a func- 
tion of 9 for each branch plotted in figure [3] shows 
a clear separation between aeolian and standard 
regimes in two different zones of 9. This is inter- 
esting since for flute-like instruments playing un- 
der normal conditions, 9 is expected to be larger 
than 4 [8]. Figure 0] shows the same trend for 
the standard regime. On the contrary, we observe 
that all aeolian regimes are located in a zone de- 
fined by 9 < 4. This feature is in agreement with 
the fact that they occur at low jet velocities. In- 
deed, for the same oscillation frequency, 9 is neces- 
sarily lower for an aeolian regime than for a stan- 
dard regime. 

Moreover, the bifurcation diagram displayed 
in figure H] shows a particular amplitude evolu- 
tion for aeolian regimes. Indeed, one notices, for 
such regimes, a bell-shaped curve, whereas for the 
"standard" regime one can observe a saturation of 
the amplitude. Thus, while the standard regime 
exists when r tends to zero (and thus when the 
jet velocity Uj tends to infinity - see equation ©), 
aeolian regimes conversely exist only for restricted 
ranges of these two parameters. 

3.2.2 Frequency evolutions 

Figure [5] presents the oscillation frequency 
function of 9, along the different periodic solution 
branches shown in figure [3j Similarly to the am- 
plitude, it reveals different evolution patterns for 
aeolian and standard regimes : 



Figure 4: Bifurcation diagram of the system ([6]) 
with same parameters value as in figure [3J Ab- 
scissa: dimensionless jet velocity 9. Ordinate: 
amplitude of the oscillating variable v(t). Sym- 
bols x represent the stable parts of the branches. 
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Figure 5: Bifurcation diagram of the system 
with same parameters as in figure [3l Abscissa: di- 
mensionless jet velocity 9. Ordinate: dimension- 
less oscillation frequency 4^ (with /i = ^ the 
resonance frequency). Symbols x represent the 
stable parts of the branches. 

• For the standard regime, one observes an im- 
portant frequency deviation just after the os- 
cillation threshold. Then the frequency tends 
toward the resonance frequency f\ = uji/2tt. 

• For aeolian regimes, one observes an impor- 
tant frequency deviation just after the oscil- 
lation threshold, followed by a second fre- 
quency deviation after a plateau around the 
resonance frequency (inflection point at f\ = 
ujx/2-k). 

3.3 Experimental illustration 

Because of the important simplification of the 
model, a quantitative comparison between 
experimental and numerical results is not pos- 
sible. However, a qualitative comparison is 
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interesting, and proves that the dynamical sys- 
tem studied in this paper reproduces some key 
features of the behaviour of flute-like instruments. 



Experimentally, aeolian sounds occur for very 
low blowing pressures (which correspond to high 
values of the delay r), that is to say when the 
musician blows very gently in the instrument. It 
is why experimental observation of such sounds is 
particularly complicated. The use of a pressure- 
controlled artificial mouth (which permits to play 
the instrument without musician, and to control 
very precisely the mouth pressure) provides 
new information about the dynamics of the 
instrument [9]. The simultaneous measurement 
of the mouth pressure (related to the jet velocity 
Uj by the stationary Bernoulli equation) and the 
acoustic pressure under the labium (see figure [T]) 
allows to study the influence of some parameters 
on the characteristics of the oscillation regimes. 

Time-frequency analysis of the sound produced 
by a recorder during both an increasing and a de- 
creasing ramp of the blowing pressure (figure [6]) 
highlights several oscillation regimes (zones A and 
C) below the threshold of the principal regime cor- 
responding to the expected note (zone B). Figure 
[7] represents, for the same data, the oscillation am- 
plitude as a function of 9. As for aeolian regimes 
on the numerical bifurcation diagram (figure H| , 
we observe that one of these regimes is located in 
a zone defined by 9 < 4, whereas the principal 
regime (zone B) appears for 9 > 4. If the second 
one (regime A2-C2 in figures [6] and [7]) appears for 
higher values of 9 (about 5 < 9 < 6.5), it nev- 
ertheless can be considered as an aeolian regime. 
Indeed, it corresponds to an oscillation on the first 
resonance mode of the pipe, whereas regimes Al- 
Cl and B corresponds to oscillations on the second 
resonance mode. Consequently, the oscillation fre- 
quency being lower, the value of 9 is increased. 

Moreover, for the two sounds which appear for 
low values of the mouth pressure, the bell-shaped 
evolution of the oscillating amplitude is clearly 
visible, recalling the characteristics of aeolian 
regimes of the model. On the contrary, the 
oscillating amplitude of the principal regime 
shows a different evolution, comparable to the 
one observed in figure 0] for the standard regime. 
As far as the frequency is concerned, large 
variations can be observed in figure particu- 
larly close to the threshold of the standard regime. 




Figure 6: Time- frequency representation (top) of 
the sound of a treble recorder played by an arti- 
ficial mouth making an increasing and decreasing 
pressure ramp (bottom). G sharp fingering (4th 
octave). 
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Figure 7: Same experimental data as in figure [6) 
reprentation of the acoustic pressure amplitude 
measured in the resonator, with respect to the di- 
mensionless jet velocity 9. Letters in legend refer 
to letters in figure [6j 
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Thereby, the very simple model studied here 
produces aeolian regimes, which present partic- 
ular features which recalls some experimental ob- 
servations on a recorder played by an artificial 
mouth. Moreover, these observations on the am- 
plitude and frequency evolution patterns for aeo- 
lian regimes can recall the results of Meissner in 
the case of a cavity resonator excited by an air jet 

4 Jumps between periodic solu- 
tion branches 

In this section, we study system ([6]) as in section[3l 
with the addition of a second acoustic mode (i.e. 
two terms in equation (j3|), m = 2). As for an open 
cylindrical resonator, the resonance frequency of 
the second mode is chosen close to twice the first 
resonance frequency [2]. 

4.1 Register change and hysteresis 
phenomenon 

The use of a transfer function including two modes 
implies the existence of two standard regimes. 
Each of these regimes results from the coupling 
between one of the two acoustic modes and the 
first hydronamic mode of the jet (n = 0). 

For a second resonance frequency slightly lower 
than twice the first one, numerical continuation 
of the corresponding periodic solution branches 
leads to the bifurcation diagram shown in figure 
[TOl - (a). As previously, it represents the oscillating 
amplitude of the variable v(t) as a function of the 
dimensionless delay f. Here again, the stability 
is adressed by examining the Floquet multipliers. 
However, there is a major difference with the case 
where only one mode of the resonator is consid- 
ered: indeed, both stable and unstable portions 
appear (symbol x corresponds to stable parts of 
the branches). 

More precisely, we can distinguish three differ- 
ent zones in figure [TU]- (a), highlighted in figure [HJ 
which is a zoom of figure [TU1- (a), for < uj\t < 2: 

• the first one is defined for uj\t > 0.7; in this 
case the first standard regime (called "first 
register" in the context of musical acoustics), 
is stable. 

• for 0.1 < uj\T < 0.7, the two standard 



regimes (i.e. the first and the second reg- 
ister) are simultaneously stable. 

• in the third zone, where oj\T < 0.1, the first 
register is unstable, and the second register 
is the only stable solution of the system. 

The existence of a range of the dimensionless 
delay f (between 0.1 and 0.7) where two peri- 
odic solutions are simultaneously stable implies 
the existence of an hysteresis phenomenon. In- 
deed, starting from a value of f where the first 
register is the only stable solution, and decreasing 
the delay f (i.e. "blowing harder"), we observe 
that the system follows the branch corresponding 
to the first register until it becomes unstable (for 
f = 0.1). At this point, the system synchronizes 
with the second register. 

Starting from this new point, and increasing the 
delay r (i.e. "blowing softer"), the system follows 
the branch corresponding to the second register. 
It is only when this second branch becomes un- 
stable (for f = 0.7) that the sytem comes back to 
the first register. 

A time-domain simulation, using a Bogacki- 
Shampine method based on a third-order Runge- 
Kutta scheme |25] allows to confirm and highlight 
this hysteresis phenomenon. In order to com- 
pare the different results, the oscillating ampli- 
tudes computed with this time-domain solver, for 
both an increasing and a decreasing ramp of the 
delay f are superimposed on the corresponding 
bifurcation diagram (figure [8}. 
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Figure 8: Zoom of figure [TUl - (a): bifurcation di- 
agram of the system §§§ , obtained with numerical 
continuation. Parameter values: m = 2 ; a = 340 
; oi = 10 ; wi = 2764; a 2 = 30 ; uj 2 = 5510; 
Qi = 100 ; Q 2 = 100. Abscissa: dimensionless 
delay lo±t. Ordinate: oscillating amplitude of the 
variable v. Symbols x represent the stable parts 
of the branches. 



7 



4.2 Register change: experimental il- 
lustration 

In flute-like instruments, it is well-known by musi- 
cians that blowing harder in the instrument even- 
tually leads to a "jump" to the note an octave 
above. Figure [9] represents the frequency evolu- 
tion of the sound produced by a treble recorder 
played by the artificial mouth, during an increas- 
ing and a decreasing ramp of the alimentation 
pressure. 
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Figure 9: Playing frequency of a treble recorder 
blown by an artificial mouth, during an increasing 
(solid line) and decreasing (dashed line) blowing 
pressure ramp, showing jumps between the two 
standard regimes. F fingering (third octave). 

These experimental results illustrate the 
register change phenomenon, corresponding 
to the "jump" from a periodic solution branch 
of frequency f\ to another of frequency /2 ~ 2 • f\ . 

An hysteresis is also observed. Recalling that 
in flute- like instruments, the blowing pressure is 
related to the convection delay r of the hydrody- 
namic perturbations along the jet, these experi- 
mental results can be compared to numerical re- 
sults presented in figure [HJ It shows again that 
the system under study reproduces classical be- 
haviour of recorder-like instruments. 

5 Influence of the ratio of the 
two resonance frequencies 

Numerical continuation allows a systematic inves- 
tigation of the parameters influence on the oscilla- 
tion regimes. In this section, we focus on the influ- 
ence of the ratio of the two resonance frequencies. 
Indeed, this parameter is well-known to instru- 
ment makers and players for its influence on the 
sound characteristics |28j . Moreover, its influence 



has been proved for other musical instruments (for 
example, in |16]). 

5.1 Stability of periodic solution 
branches 

Figures [TOl present bifurcation diagrams obtained 
for different values of this ratio (respectively ^ = 
1.99 ; = 2 et ^ = 2.05). Only the second 
resonance frequency LO2 varies from one case to the 
other; all the others parameters are kept constant. 

For all these bifurcation diagrams, the bell- 
shaped curve corresponds to an aeolian regime 
(see section I3.2p , the dark gray curve corresponds 
to the first register, the light gray dashed one cor- 
responds to the second register, and symbols x 
represent stable parts of the branches. 

Although the resonance frequency of the sec- 
ond mode is only slightly modified from one case 
to the other, one observes significant changes in 
the stability properties of the two registers. In 
the case of two perfectly harmonic acoustic modes 
(u)2 = 2wi), represented in figure [TOl- (b), the first 
register is initially stable, and becomes unstable 
when the delay f decreases (i.e. when the blowing 
pressure increases). On the contrary, the second 
register is first unstable and becomes stable when 
the delay f decreases. For the range of f located 
between the loss of stability of the first register 
and the stabilization of the second register (i.e. 
between the two vertical dot-dashed lines in fig- 
ure [10]- (b)), there is no stable periodic or static 
solution. 

In order to determine the nature of the differ- 
ent resulting bifurcations, we propose to study 
Floquet multipliers in the vicinity of stability 
changes. As shown in figure [UJ stabilization of 
the second register occurs, at f = 0.68, by the 
crossing of the unit circle by two conjugate com- 
plex Floquet multipliers. 

In the case of a direct bifurcation, this Hopf 
bifurcation of the limit cycle leads to a quasi pe- 
riodic regime |10j . The numerical tools used here 
do not allow numerical continuation of this kind 
of solutions. On the other hand, time-domain 
integration methods allow to study such regimes. 
Figure [12] shows the time-frequency analysis of 
the signal obtained with the same time-domain 
solver as previously, by increasing the dimension- 
less delay f in a quasi-static way from an initial 
value for which the second register is stable, 
to a final value for which the first register is 
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climensionle^s delay x ro (w = 276tf" 



500p 

> 

Q) 

(0 
> 

I 300 
•a 

3 

I 200 " 
O) 

c 

I 100- 



static solutions 
— standard regime (1 st register, n = 0) 

standard regime (2 nd register, n = 0) 
— aeolian regime (n = 1 ) 
x stable solutions 




dimensionless delay t to (ra = 2764f 



u 

o- -k 










static solutions 

standard regime (1 st register, n = 0) 

„_■„. ,„nd ._ 










sianaara regime \t register, n = u) 
— aeolian regime (n = 1) 
x stable solutions 


o- \ 




























o- 














o 1 : 


3 L 
imensionless delay x a 




f 


7 



Figure 10: Bifurcation diagrams of the system ([6]), 
obtained with numerical continuation, for differ- 
ent values of inharmonicity of the two resonance 
modes. Parameter values: m = 2 ; a = 340 ; 
ax = 10 ; a 2 = 30 ; Q 1 = 100 ; Q 2 = 100 ; 
wi = 2764. (a) Top: lo 2 = 5510, ^ = 1.99. 
(b) Middle: u 2 = 5528, ^ = 2. (c) Bottom: 
oj 2 = 5654, ^ = 2.05. Abscissa: dimensionless 
delay uj\t. Ordinate: oscillating amplitude of the 
variable v. Symbols x represent the stable parts 
of the branches. 



stable. We note three different regimes: zones 
A {uj\t < 0.71) and C {uj\t > 0.9) correspond 
respectively to the second and first registers. 
In zone B (0.71<WiT < 0.9), the presence of 
multiple frequencies reveals the existence of a 
quasiperiodic regime, which agrees with the 
results of Floquet analysis. 
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Figure 11: Representation in the complex plane 
of the Floquet multipliers for a point located just 
before the stabilization of the second periodic so- 
lution branch (related to the second register) . Pa- 
rameter values are the same as for figure [TO] - (b) . 



7000 



6000 



5000 



t4000 



3 3000 



2000 



1000 



3.6 0.65 0.7 _,. 0.75 , 0,8 , 0.85 „ 0.9 0.95 
dimensionless delay wi ("•= 2764) 



-5i 



-50 
-10( 
I|l5( 
|-20( 



Figure 12: Time-frequency representation ob- 
tained by a time-domain simulation of the system 
©, using a third-order Runge-Kutta scheme. A 
slowly increasing ramp of the dimensionless delay 
uj\T is achieved. Zone A corresponds to the second 
register, zone C corresponds to the first register, 
and zone B to a quasiperiodic regime. Parameter 
values are the same as for figure [TO] - (b) . 



As shown in figure [TO] - (a), a small decrease 
of the second mode resonance frequency (of 
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about 0.6 %) compared to the harmonic case 
0J2 = 2wi does not alter the general form of the 
bifurcation diagram: the first register is stable 
for high values, whereas the second register is 
first unstable and becomes stable when the delay 
f decreases. However, the stability ranges are 
significantly modified. Compared to the previous 
case, the first register is stable for smaller values 
of the delay f, which leads to the existence of a 
range of f where the two registers are simulta- 
neously stable. In the later case, the hysteresis 
phenomenon emphasized in section [5] becomes 
possible. 

These stability ranges are affected again by an 
increase of the second mode resonance frequency 
(leading to a difference of about 2% between 0J2 
and 2wi). Figure E3 - (c) shows the resulting 
bifurcation diagram. One can notice that the first 
register is now stable for all values of f between 
2 and 0. Hence, once the system is synchronized 
on this branch of solutions, no register change is 
possible, and the only way to make the system 
oscillate on the second register is to change initial 
conditions. Time-domain simulations results (not 
presented here) show good agreement with these 
characteristics. 



5.2 Experimental illustration 

As emphasized in section E the simplicity of the 
studied system prevents a quantitative compari- 
son between numerical results and experimental 
datas. However, we can qualitatively compare 
some experimental phenomena with numerical 
results presented in the previous section. 

The use of an artificial mouth with a real instru- 
ment highlights the influence of the inharmonicity 
on the sound characteristics. Indeed, the mea- 
sured inharmonicity depends on the fingering. 

Thus, for a small inharmonicity ~ 2.04), an 
increasing mouth pressure ramp (corresponding 
to a decrease of the convection delay r) causes a 
regime change from the first register to the second 
register, including an hysteresis effect (as shown 
in figure [TO]) . 

For a larger inharmonicity of the resonance 
frequencies (^ ~ 1.92), the same experiment 
leads to a very different behaviour. As shown 
in figure Q3J we observe a transition from the 




time (s) 

Figure 13: Time-frequency representation (top) 
of the sound of a treble recorder played by an ar- 
tificial mouth during an increasing and decreas- 
ing blowing pressure ramp (bottom). G fingering 
(3th octave), with a slightly positive inharmonic- 
ity (^-2.04). 

first register to a quasiperiodic regime whitout 
hysteresis effect. 




time (s) 

Figure 14: Time-frequency representation (top) of 
the sound of a treble recorder played by an arti- 
ficial mouth during an increasing and decreasing 
blowing pressure ramp (bottom). B-flat finger- 
ing (3th octave), with a significant negative in- 
harmonicity (^ « 1.92). 

Comparison of these results with those of fig- 
ures [TO] - (a), [TOJ- (b) and [TOJ- (c) proves that the 
behaviour of both the studied system and the real 
instrument strongly depends on the inharmonicity 
of the two first resonance frequencies. 

A small change of this parameter value is 
enough to alter the oscillation thresholds of 
the different regimes, their stability properties, 
or even the nature of the oscillation regimes. 
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Depending on the case, the system can "jump" 
between branches of periodic solutions (with hys- 
teresis effect), or jump from a periodic solution 
branch to a quasiperiodic regime. 

Moreover, these results seem to be consistent 
with the experience of a recorder maker [28j: ac- 
cording to his experience, the first register is sta- 
ble for a wider range of alimentation pressure in 
the case of strong inharmonicity and the case of a 
perfect harmonicity should be avoided to prevent 
instabilities. 

6 Conclusion 

We have studied a nonlinear delay dynamical 
sytem with small number of degrees of freedom. 
We focused on periodic solutions. Indeed, this 
system is inspired by the Sanctioning of flute-like 
instruments and in most cases, notes produced 
by these instruments are periodic oscillations. 

Because of the drastic simplifications, the 
studied system can not be considered a priori as 
a model of the functioning of these instruments. 
However, as shown in this paper, the use of 
numerical continuation tools, providing a more 
global vision of the dynamics of the system, 
highlights that the system not only presents a 
considerable variety of periodic regimes, but also 
has a second interest. Indeed, it qualitatively 
reproduces a lot of phenomena observed exper- 
imentally on flute-like instruments : amplitude 
and frequency evolutions for both standard and 
aeolian regimes, regime changes with hysteresis 
effect, and quasi periodic oscillations. 

We can furthermore investigate the influence of 
some parameters related to instrument maker's 
issues. We focused here on the role of the inhar- 
monicity of the two first resonance frequencies. 
Analysis of bifurcation diagrams leads to results 
that are consistent with the empirical knowledge 
of an instrument maker. 

However, it would be hazardous to use the stud- 
ied dynamical system predictive tool, for ex- 
ample for the design of musical instruments. In- 
deed, some parameters of the system can not be 
related to mesurable physical quantities. More- 
over, as we emphasized in section [TJ some impor- 
tant physical phenomena have not been taken into 



account. These two elements make vain any at- 
tempt of quantitative comparison between numer- 
ical and experimental results. Particularly, we did 
not discussed about the sound timbre, which cor- 
responds to the spectrum of the sound produced 
by the instrument. This characteristics is essen- 
tial in a musical context since it allows us to dis- 
tinguish, in the case of a steady-state regime, the 
sound of a flute from that, for example, of a trum- 
pet. The system studied here can not be used 
to predict the sound timbre, since some elements 
known for their significant influence on the spec- 
trum have been neglected: 

• the offset between the position of the labium 
and the channel exit. Due to symmetry prop- 
erties of the nonlinear function, this parame- 
ter controls the ratio between even and odd 
harmonics. [26] 

• nonlinear losses due to the flow separation at 
the labium, which have an important influ- 
ence on the sound level. [23] 

• the dipolar character of the pressure source, 
created by the oscillation of the jet around 
the labium, which favors high frequencies 

mi- 

Taking into account these two first elements do 
not compromise the use of the approach presented 
in this paper. Only the number of parameters 
and degree of freedom (and thus the computation 
time) increases. However, taking into account the 
third element involves the presence of a delayed 
derivative term in the right-hand side of the sec- 
ond equation of system ([JJ). This change trans- 
forms the delay system into a neutral delay dy- 
namical system, and thus introduces additional 
difficulties in numerical continuation and system 
analysis |17| . 
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A DDE-Biftool: Reformulation 
of the studied system 

To be implemented in the software DDE-Biftool, 
system (JTJ) needs to be reformulated as follows: 



x = f{x(t),x{t - r),7) 



(9) 



where r is a delay, x the vector of state variables, 
and 7 the set of parameters. 



Equation (|15p becomes: 



It leads to: 



(16) 



We detail here such a reformulation in the case of 
a single- mode transfer function (equation jl])): 



juj ■ ai 



i )2 , ,2 i „■, , m 



(10) 



where u) is the pulsation, and a±, uj\ and Qi are 
respectively the modal amplitude, the resonance 
pulsation and the quality factor of the resonance 
mode. 



Let us inject equation ()10p in the third equation 
of system ([1]): 



V ■ 



VI 



juj ■ ai ■ P. 



An inverse Fourier transform leads to: 

v(t) + uj ■ v{t) + ^- • v(t) = at • p(i). (12) 
VI 

The right-hand side of this expression can be de- 
veloped using the second equation of system ([T]): 



p(t) = a ■ tanh [z(t)\ 



(13) 



Ht) + -^r-v(t)+v(t) 



}{t-f) ■ {l - tanh 2 [v{t-f)}} , 



(17) 



where v define the derivative of v with respect to 
dimensionless time t. 
We define a new variable: 



y(t) = v(t), 



(18) 



which leads to the correct form of the system 
(given by equation d!])): 

{ v(t)=y(t) 

~ = aia Q_ f y^_ tanh 2 [ v (t-?)] } 

"(t) - 7^-l/(t)- 

(19) 



Generalization to a system including m res- 
onance modes is straightforward. In this case, 
equation contains m terms (with m > 1): 



Differentiating with respect to time, and using the 
first equation of system ([1]), we obtain: 



p(t) = a ■ z(t) • {1 - tanh 2 [z(t)}} 



(14) 



;(i -r) ■ {1 - tanh 2 [v(t - r)]} 
Injecting this expression in equation (jl2[) leads to: 



a • in 



i)(i) + -pr-v(t) + ui 2 v(t) = a\a ■ v(t — r) 
VI 

• {l - tanh 2 [v(t - r)]} 
(15) 

To improve numerical conditioning of the prob- 
lem, it is convenient to make the temporal vari- 
ables dimensionless. Let us introduce the new 
variables: 

i = oj\ • t 

T = U\-T 



Yin — ^2 



ju ■ a k 



k=l k J Qk 



(20) 



and it leads to a system of 2m equations (i.e. 2 
equations by resonance mode), of the following 
form: 



Vi = [1, 2..., m] : 

.~. OjQ I n 

yi\t) = \ 1 ~~ tanh 



,fc=l 



m / \ 2 

EN*- f )]- - *® 



■I/i(t). 



(21) 
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Figure 15: Real part of the eigenvalues A of the 
Jacobian matrix of system ([6]) along the static so- 
lution branch, with respect to the dimensionless 
delay r • uj\ . Hopf bifurcations (intersections with 
the horizontal line Re{\) = 0) are highlighted 
by circles. Parameters value: m = 1, a = 10 ; 
ax = 70 ; uji = 2260 ; Qi = 50. 

B Linear analysis around the 
equilibrium solution 

B.l Eigenvalues of the Jacobian 

Linearisation of system ([6]) around equilibrium so- 
lution ([7|) allows to compute the eigenvalues A of 
the Jacobian matrix of the studied system along 
this branch. Their real parts Re{\) determine the 
stability properties of the static solution: it is sta- 
ble if and only if all the eigenvalues have negative 
real parts |10| . 

Considering, for sake of clarity, a single acous- 
tic mode of the resonator (i.e. m = 1 in equation 
dU), we represent, in figure [151 -Re(A) with re- 
spect to the continuation parameter f, along the 
branch of equilibrium solutions. This represen- 
tation highlights that this static solution is sta- 
ble for 1.8 < f < 4.1 and for 9.1 < f < 9.5. 
Analysis of both real and imaginary parts of the 
eigenvalues shows that the equilibrium solution 
is a focus point when it is stable, and a saddle 
point when it is unstable |10j . Moreover, inspec- 
tion of the intersection points with the x-axis (de- 
fined by Re(X) = 0) determines Hopf bifurcations 
jlOj . which correspond to the birth of the differ- 
ent periodic solutions branches (here highlighted 
by circles at f = 1.8 ; f = 4.1 ; f = 9.1 and 
f = 9.5). 

B.2 Open- loop gain 

In a feedback loop system such as the one pre- 
sented in section [2j another approach to get in- 



formation about the destabilization mechanisms 
of the equilibrium solution is to study the open- 
loop gain G i(u) of the linearised system. Indeed, 
the emergence of auto-oscillations from the equi- 
librium solution (corresponding to a Hopf bifur- 
cation) is possible under the two following condi- 
tions |15| : 

• modulus G of the linearised open-loop gain 
G i must be equal to 1. 

• phase P of the linearised open-loop gain G Q i 
must be a multiple of 2tt. 

Writting system ([1]) in the frequency domain 
leads to: 

Z(u) = V(co)e~ juJT 
P(u) = a ■ tanh{Z{uj)) (22) 
V{J) = Y(u) ■ P{oj) 

and linearisation of the second equation around 
the equilibrium solution ([7]) leads to the open-loop 
gain: 

G ol (u) = aY (cj)-e^, (23) 

The conditions of emergence of auto-oscillations 
are then given by: 

f G(u) = a-\Y(u J )\ = l 

\ P{uj)=arg{Y{u))-UT = -n-2Ti K ' 

where n is an integer. 

To exemplify the conditions given by equation 
()24p . let us consider a transfer function Y[u) 
representing the first two modes of a cylindri- 
cal resonator (m = 2 in equation ([J}). Figure 
[TBI shows the variables G (top) and P (bottom) 
with respect to u. Two different values of the 
delay r are considered (corresponding to two dif- 
ferent values of the jet velocity Uj = 6.5ms -1 and 
Uj = 15.7ms" 1 ) for phase P (G is independent 
of t). The points where equations (|24p are ful- 
filled are marked with circles (o) when n = and 
squares (□) when n = 1. Plots of P for lower val- 
ues of Uj (i.e. for higher values of t) would have 
revealed other intersections with horizontal lines 
P = —n ■ 2tt, with larger values of n. 

Provided the amplification a is large enough, 
this example highlights the existence of an infinity 
of solutions of equations (|24p for different values 
of Uj, each solution being related to a given value 
of the integer n. Therefore, for each value of n, 
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Figure 16: Open-loop gain modulus G (top) and 
phase P (bottom) defined by eq. (I24p , for two dif- 
ferent values of the jet velocity Uj. Emergence of 
self-sustained oscillations is possible if the mod- 
ulus is equal to one (dash line, at the top), and 
if the phase crosses a straight line with equation 
P = —n ■ 2tt (dash lines, at the bottom). 

an instability may emerge from the equilibrium 
branch. 

From the physics point of view, the integer n 
represents the rank of the hydrodynamic mode of 
the jet involved in the emergence of the instability. 
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